####################################################################
#***2019 analysis
## Packages
# To install and open the R packages that you need for this code. 
need <- c('tidyverse','readstata13','lfe','glue','rdrobust', 'stargazer','arm', 'broom', 'ggplot2', 'dotwhisker', 'gridExtra')
have <- need %in% rownames(installed.packages()) 
if(any(!have)) install.packages(need[!have]) 
invisible(lapply(need, library, character.only=T)) 

# Change path to whereever you place the models
# To set up the working directory. 
script_folder = dirname(rstudioapi::getSourceEditorContext()$path)
setwd(glue('{script_folder}'))
rm(list = ls())
setwd("../")

## load in RD data:
load("5prepdata/final_rd_data.RData")
load("5prepdata/final_rd_data_with_opp.RData")
# estimates by nmber of terms running

main.specification <- function(x) {
  rd.data$cutoff<-NULL
  rd.data$bandwidth<-NULL
  rd.data$kw<-NULL
  
  y<-rdbwselect(x, rd.data$victory_marg, p=1, c = 0,  kernel = "tri", bwselect="mserd", covs = rd.data$running_terms_served + rd.data$dem + rd.data$rep + rd.data$spc_elec + rd.data$term_length + rd.data$number_candidates + rd.data$year)
  rd.data$bandwidth<-y$bws[1]
  rm(y)
  
  # gen cutoff = 0
  rd.data$cutoff <- 0
  # gen kw = 1-(abs(cutoff-victory_marg))/bandwidth
  rd.data$kw <- 1-(abs(rd.data$cutoff-rd.data$victory_marg))/rd.data$bandwidth
  # replace kw = 0 if (victory_marg>cutoff+bandwidth | victory_marg<cutoff-bandwidth) 
  rd.data$kw[rd.data$victory_marg>(rd.data$cutoff+rd.data$bandwidth) | rd.data$victory_marg<(rd.data$cutoff-rd.data$bandwidth)] <- 0
  
  out <- felm(x[rd.data$kw>0] ~ won_election + victory_marg + f_x_win + running_terms_served + dem + rep + spc_elec + term_length + number_candidates | year + state | 0 | state, data = rd.data[rd.data$kw>0,], weights = rd.data$kw[rd.data$kw>0])
  return(out)
}

main.specification.opp <- function(x) {
  data.orig$cutoff<-NULL
  data.orig$bandwidth<-NULL
  data.orig$kw<-NULL
  
  y<-rdbwselect(x, data.orig$victory_marg, p=1, c = 0,  kernel = "tri", bwselect="mserd", covs = data.orig$running_terms_served + data.orig$dem + data.orig$rep + data.orig$spc_elec + data.orig$term_length + data.orig$number_candidates + data.orig$year)
  data.orig$bandwidth<-y$bws[1]
  rm(y)
  
  # gen cutoff = 0
  data.orig$cutoff <- 0
  # gen kw = 1-(abs(cutoff-victory_marg))/bandwidth
  data.orig$kw <- 1-(abs(data.orig$cutoff-data.orig$victory_marg))/data.orig$bandwidth
  # replace kw = 0 if (victory_marg>cutoff+bandwidth | victory_marg<cutoff-bandwidth) 
  data.orig$kw[data.orig$victory_marg>(data.orig$cutoff+data.orig$bandwidth) | data.orig$victory_marg<(data.orig$cutoff-data.orig$bandwidth)] <- 0
  
  out <- felm(x[data.orig$kw>0] ~ won_election + victory_marg + f_x_win + running_terms_served + dem + rep + spc_elec + term_length + number_candidates | year + state | 0 | state, data = data.orig[data.orig$kw>0,], weights = data.orig$kw[data.orig$kw>0])
  return(out)
}


## flag runhouse
# create variables first
data.orig$flag_1_runhouse <- data.orig$opp_1_runhouse
data.orig$flag_2_runhouse <- data.orig$opp_2_runhouse - data.orig$opp_1_runhouse
data.orig$flag_3_runhouse <- data.orig$opp_3_runhouse - data.orig$opp_2_runhouse
data.orig$flag_4_runhouse <- data.orig$opp_4_runhouse - data.orig$opp_3_runhouse
data.orig$flag_5_runhouse <- data.orig$opp_5_runhouse - data.orig$opp_4_runhouse
data.orig$flag_6_runhouse <- data.orig$opp_6_runhouse - data.orig$opp_5_runhouse
data.orig$flag_7_runhouse <- data.orig$opp_7_runhouse - data.orig$opp_6_runhouse
data.orig$flag_8_runhouse <- data.orig$opp_8_runhouse - data.orig$opp_7_runhouse
data.orig$flag_9_runhouse <- data.orig$opp_9_runhouse - data.orig$opp_8_runhouse
data.orig$flag_10_runhouse <- data.orig$opp_10_runhouse - data.orig$opp_9_runhouse
data.orig$flag_11_runhouse <- data.orig$opp_11_runhouse - data.orig$opp_10_runhouse
data.orig$flag_12_runhouse <- data.orig$opp_12_runhouse - data.orig$opp_11_runhouse
data.orig$flag_13_runhouse <- data.orig$opp_13_runhouse - data.orig$opp_12_runhouse
data.orig$flag_14_runhouse <- data.orig$opp_14_runhouse - data.orig$opp_13_runhouse
data.orig$flag_15_runhouse <- data.orig$opp_15_runhouse - data.orig$opp_14_runhouse



#flag_1_hle etc
m1 <- main.specification(data.orig$flag_1_runhouse) %>% tidy %>% filter(term == "won_election")
m2 <- main.specification(data.orig$flag_2_runhouse) %>% tidy %>% filter(term == "won_election")
m3 <- main.specification(data.orig$flag_3_runhouse) %>% tidy %>% filter(term == "won_election")
m4 <- main.specification(data.orig$flag_4_runhouse) %>% tidy %>% filter(term == "won_election")
m5 <- main.specification(data.orig$flag_5_runhouse) %>% tidy %>% filter(term == "won_election")
m6 <- main.specification(data.orig$flag_6_runhouse) %>% tidy %>% filter(term == "won_election")
m7 <- main.specification(data.orig$flag_7_runhouse) %>% tidy %>% filter(term == "won_election")
m8 <- main.specification(data.orig$flag_8_runhouse) %>% tidy %>% filter(term == "won_election")
m9 <- main.specification(data.orig$flag_9_runhouse) %>% tidy %>% filter(term == "won_election")
m10 <- main.specification(data.orig$flag_10_runhouse) %>% tidy %>% filter(term == "won_election")
m11 <- main.specification(data.orig$flag_11_runhouse) %>% tidy %>% filter(term == "won_election")
m12 <- main.specification(data.orig$flag_12_runhouse) %>% tidy %>% filter(term == "won_election")
m13 <- main.specification(data.orig$flag_13_runhouse) %>% tidy %>% filter(term == "won_election")
m14 <- main.specification(data.orig$flag_14_runhouse) %>% tidy %>% filter(term == "won_election")
m15 <- main.specification(data.orig$flag_15_runhouse) %>% tidy %>% filter(term == "won_election")

m1$model <- "1"
m2$model <- "2"
m3$model <- "3"
m4$model <- "4"
m5$model <- "5"
m6$model <- "6"
m7$model <- "7"
m8$model <- "8"
m9$model <- "9"
m10$model <- "10"
m11$model <- "11"
m12$model <- "12"
m13$model <- "13"
m14$model <- "14"
m15$model <- "15"




m.hle <- bind_rows(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15)
#dwplot(m.hle)

m.hle$model <- as.numeric(m.hle$model)

m.hle$ub <- m.hle$estimate + 2*m.hle$std.error
m.hle$lb <- m.hle$estimate -2* m.hle$std.error

ggplot(m.hle, aes(x = model, y = estimate)) + geom_point() +
  geom_errorbar(aes(ymin = lb, ymax=ub), width=0.2) +
  scale_x_continuous(breaks = c(seq(1, 15)), labels = c(seq(1, 15))) + 
  #ylim(-0.005, 0.025) +
  geom_hline(yintercept = 0, linetype=2, color="#666666") +
  geom_line() +
  theme_bw() +
  theme(panel.grid.minor.y=element_blank(),
        panel.grid.minor.x = element_blank()) +
  xlab("Specific opportunity window") + ylab("Estimate") 
  ggsave(file = "7tex/manuscript/tables/sourcefiles/Figure 5c.pdf", units = "in", height=4, width=6)


### ever run house
m1 <- main.specification(rd.data$flag_1_winhouse) %>% tidy %>% filter(term == "won_election")
m2 <- main.specification(rd.data$flag_2_winhouse) %>% tidy %>% filter(term == "won_election")
m3 <- main.specification(rd.data$flag_3_winhouse) %>% tidy %>% filter(term == "won_election")
m4 <- main.specification(rd.data$flag_4_winhouse) %>% tidy %>% filter(term == "won_election")
m5 <- main.specification(rd.data$flag_5_winhouse) %>% tidy %>% filter(term == "won_election")
m6 <- main.specification(rd.data$flag_6_winhouse) %>% tidy %>% filter(term == "won_election")
m7 <- main.specification(rd.data$flag_7_winhouse) %>% tidy %>% filter(term == "won_election")
m8 <- main.specification(rd.data$flag_8_winhouse) %>% tidy %>% filter(term == "won_election")
m9 <- main.specification(rd.data$flag_9_winhouse) %>% tidy %>% filter(term == "won_election")
m10 <- main.specification(rd.data$flag_10_winhouse) %>% tidy %>% filter(term == "won_election")
m11 <- main.specification(rd.data$flag_11_winhouse) %>% tidy %>% filter(term == "won_election")
m12 <- main.specification(rd.data$flag_12_winhouse) %>% tidy %>% filter(term == "won_election")
m13 <- main.specification(rd.data$flag_13_winhouse) %>% tidy %>% filter(term == "won_election")
m14 <- main.specification(rd.data$flag_14_winhouse) %>% tidy %>% filter(term == "won_election")
m15 <- main.specification(rd.data$flag_15_winhouse) %>% tidy %>% filter(term == "won_election")

m1$model <- "1"
m2$model <- "2"
m3$model <- "3"
m4$model <- "4"
m5$model <- "5"
m6$model <- "6"
m7$model <- "7"
m8$model <- "8"
m9$model <- "9"
m10$model <- "10"
m11$model <- "11"
m12$model <- "12"
m13$model <- "13"
m14$model <- "14"
m15$model <- "15"


m.winhouse <- bind_rows(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15)
#dwplot(m.hle)

m.winhouse$model <- as.numeric(m.winhouse$model)

m.winhouse$ub <- m.winhouse$estimate + 2*m.winhouse$std.error
m.winhouse$lb <- m.winhouse$estimate -2* m.winhouse$std.error

ggplot(m.winhouse, aes(x = model, y = estimate)) + geom_point() +
  geom_errorbar(aes(ymin = lb, ymax=ub), width=0.2) +
  scale_x_continuous(breaks = c(seq(1, 15)), labels = c(seq(1, 15))) + 
  #ylim(-0.005, 0.025) +
  geom_hline(yintercept = 0, linetype=2, color="#666666") +
  geom_line() +
  theme_bw() +
  theme(panel.grid.minor.y=element_blank(),
        panel.grid.minor.x = element_blank()) +
  xlab("Specific opportunity window") + ylab("Estimate") 
  ggsave(file = "7tex/manuscript/tables/sourcefiles/Figure 5d.pdf", units="in", height=4, width=6)



#### opportunity to run windows ####

m1 <- main.specification.opp(data.orig$opp_1_runhouse) %>% tidy %>% filter(term == "won_election")
m2 <- main.specification.opp(data.orig$opp_2_runhouse) %>% tidy %>% filter(term == "won_election")
m3 <- main.specification.opp(data.orig$opp_3_runhouse) %>% tidy %>% filter(term == "won_election")
m4 <- main.specification.opp(data.orig$opp_4_runhouse) %>% tidy %>% filter(term == "won_election")
m5 <- main.specification.opp(data.orig$opp_5_runhouse) %>% tidy %>% filter(term == "won_election")
m6 <- main.specification.opp(data.orig$opp_6_runhouse) %>% tidy %>% filter(term == "won_election")
m7 <- main.specification.opp(data.orig$opp_7_runhouse) %>% tidy %>% filter(term == "won_election")
m8 <- main.specification.opp(data.orig$opp_8_runhouse) %>% tidy %>% filter(term == "won_election")
m9 <- main.specification.opp(data.orig$opp_9_runhouse) %>% tidy %>% filter(term == "won_election")
m10 <- main.specification.opp(data.orig$opp_10_runhouse) %>% tidy %>% filter(term == "won_election")
m11 <- main.specification.opp(data.orig$opp_11_runhouse) %>% tidy %>% filter(term == "won_election")
m12 <- main.specification.opp(data.orig$opp_12_runhouse) %>% tidy %>% filter(term == "won_election")
m13 <- main.specification.opp(data.orig$opp_13_runhouse) %>% tidy %>% filter(term == "won_election")
m14 <- main.specification.opp(data.orig$opp_14_runhouse) %>% tidy %>% filter(term == "won_election")
m15 <- main.specification.opp(data.orig$opp_15_runhouse) %>% tidy %>% filter(term == "won_election")



m1$model <- "1"
m2$model <- "2"
m3$model <- "3"
m4$model <- "4"
m5$model <- "5"
m6$model <- "6"
m7$model <- "7"
m8$model <- "8"
m9$model <- "9"
m10$model <- "10"
m11$model <- "11"
m12$model <- "12"
m13$model <- "13"
m14$model <- "14"
m15$model <- "15"




m.hle <- bind_rows(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15)
#dwplot(m.hle)

m.hle$model <- as.numeric(m.hle$model)

m.hle$ub <- m.hle$estimate + 2*m.hle$std.error
m.hle$lb <- m.hle$estimate -2* m.hle$std.error

ggplot(m.hle, aes(x = model, y = estimate)) + geom_point() +
  geom_errorbar(aes(ymin = lb, ymax=ub), width=0.2) +
  scale_x_continuous(breaks = c(seq(1, 15)), labels = c(seq(1, 14), "Ever")) + 
  #ylim(-0.005, 0.025) +
  geom_hline(yintercept = 0, linetype=2, color="#666666") +
  geom_line() +
  theme_bw() +
  theme(panel.grid.minor.y=element_blank(),
        panel.grid.minor.x = element_blank()) +
  xlab("Culumative opportunities") + ylab("Estimate") 
  ggsave(file = "7tex/manuscript/tables/sourcefiles/Figure 5a.pdf", units = "in", height=4, width=6)


## win house

m1 <- main.specification.opp(data.orig$opp_1_winhouse) %>% tidy %>% filter(term == "won_election")
m2 <- main.specification.opp(data.orig$opp_2_winhouse) %>% tidy %>% filter(term == "won_election")
m3 <- main.specification.opp(data.orig$opp_3_winhouse) %>% tidy %>% filter(term == "won_election")
m4 <- main.specification.opp(data.orig$opp_4_winhouse) %>% tidy %>% filter(term == "won_election")
m5 <- main.specification.opp(data.orig$opp_5_winhouse) %>% tidy %>% filter(term == "won_election")
m6 <- main.specification.opp(data.orig$opp_6_winhouse) %>% tidy %>% filter(term == "won_election")
m7 <- main.specification.opp(data.orig$opp_7_winhouse) %>% tidy %>% filter(term == "won_election")
m8 <- main.specification.opp(data.orig$opp_8_winhouse) %>% tidy %>% filter(term == "won_election")
m9 <- main.specification.opp(data.orig$opp_9_winhouse) %>% tidy %>% filter(term == "won_election")
m10 <- main.specification.opp(data.orig$opp_10_winhouse) %>% tidy %>% filter(term == "won_election")
m11 <- main.specification.opp(data.orig$opp_11_winhouse) %>% tidy %>% filter(term == "won_election")
m12 <- main.specification.opp(data.orig$opp_12_winhouse) %>% tidy %>% filter(term == "won_election")
m13 <- main.specification.opp(data.orig$opp_13_winhouse) %>% tidy %>% filter(term == "won_election")
m14 <- main.specification.opp(data.orig$opp_14_winhouse) %>% tidy %>% filter(term == "won_election")
m15 <- main.specification.opp(data.orig$opp_15_winhouse) %>% tidy %>% filter(term == "won_election")



m1$model <- "1"
m2$model <- "2"
m3$model <- "3"
m4$model <- "4"
m5$model <- "5"
m6$model <- "6"
m7$model <- "7"
m8$model <- "8"
m9$model <- "9"
m10$model <- "10"
m11$model <- "11"
m12$model <- "12"
m13$model <- "13"
m14$model <- "14"
m15$model <- "15"




m.hle <- bind_rows(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15)
#dwplot(m.hle)

m.hle$model <- as.numeric(m.hle$model)

m.hle$ub <- m.hle$estimate + 2*m.hle$std.error
m.hle$lb <- m.hle$estimate -2* m.hle$std.error

ggplot(m.hle, aes(x = model, y = estimate)) + geom_point() +
  geom_errorbar(aes(ymin = lb, ymax=ub), width=0.2) +
  scale_x_continuous(breaks = c(seq(1, 15)), labels = c(seq(1, 14), "Ever")) + 
  #ylim(-0.005, 0.025) +
  geom_hline(yintercept = 0, linetype=2, color="#666666") +
  geom_line() +
  theme_bw() +
  theme(panel.grid.minor.y=element_blank(),
        panel.grid.minor.x = element_blank()) +
  xlab("Cululative opportunities") + ylab("Estimate") 
  ggsave(file = "7tex/manuscript/tables/sourcefiles/Figure 5b.pdf", units = "in", height=4, width=6)













#### empirical RD estimates ####


##  THIS FUNCTION PRODUCES FIGURES 3 AND 4A AND 4B, AS WELL AS APPENDIX FIGURES 2A-D
rd.plot <- function(band.var, yvar){
  y<-rdbwselect(band.var, rd.data$victory_marg, p=1, c = 0,  kernel = "tri", bwselect="mserd", covs = rd.data$running_terms_served + rd.data$dem + rd.data$rep + rd.data$spc_elec + rd.data$term_length + rd.data$number_candidates + rd.data$year)
  rd.data$bandwidth<-y$bws[1]
  
  file.name <- paste("7tex/manuscript/tables/sourcefiles/rdplot_", quo_name(enquo(yvar)), ".pdf", sep="")
 
  yvar <- enquo(yvar)
  
  point.data <- rd.data %>%
    filter(rd.data$victory_marg > -rd.data$bandwidth[1] & rd.data$victory_marg < rd.data$bandwidth[1]) %>%
    mutate(round.marg = round(victory_marg, digits = 2)) %>%
    group_by(round.marg) %>%
    summarize(hle.prob = mean(!! yvar))
  
    rd.data %>% ungroup %>%
    filter(rd.data$victory_marg > -rd.data$bandwidth[1] & rd.data$victory_marg < rd.data$bandwidth[1]) %>%
    mutate(win = ifelse(victory_marg > 0, "win", "lose")) %>%
    ggplot(aes(x = victory_marg, y = !!yvar, group = win)) +
    geom_smooth(method="lm",  formula=y ~ poly(x, 2)) +
    geom_point(data = point.data, aes(x = round.marg, y = hle.prob), inherit.aes = F, alpha=.5) +
    geom_vline(xintercept = 0, linetype = 2, color = "#666666", alpha=.8) +
    theme_bw() +
    theme(panel.grid.minor.y=element_blank(),
          panel.grid.minor.x = element_blank()) +
    xlab("Victory Margin") + ylab("Probability") 
    #ggsave(file = file.name, units="in", width=6, height=4)
}


rd.plot(rd.data$ever_hle, ever_hle)

rd.plot(rd.data$ever_ranprim_h, ever_ranprim_h) 
ggsave("7tex/manuscript/tables/sourcefiles/Figure2a.pdf", width=6, height=4)

rd.plot(rd.data$ever_ranprim_s, ever_ranprim_s)
ggsave("7tex/manuscript/tables/sourcefiles/Appendix Figure 2a.pdf", width=6, height=4)

rd.plot(rd.data$ever_runhouse, ever_runhouse)
ggsave("7tex/manuscript/tables/sourcefiles/Figure3.pdf", width=6, height=4)


rd.plot(rd.data$ever_runsenate, ever_runsenate)
ggsave("7tex/manuscript/tables/sourcefiles/Appendix Figure 2c.pdf", width=6, height=4)

rd.plot(rd.data$ever_gen, ever_gen)

rd.plot(rd.data$ever_winhle, ever_winhle)

rd.plot(rd.data$ever_wingen, ever_wingen)
ggsave("7tex/manuscript/tables/sourcefiles/Figure4b.pdf", width=6, height=4)

rd.plot(rd.data$ever_winprim_h, ever_winprim_h)
ggsave("7tex/manuscript/tables/sourcefiles/Figure2b.pdf", width=6, height=4)


rd.plot(rd.data$ever_winprim_s, ever_winprim_s)
ggsave("7tex/manuscript/tables/sourcefiles/Appendix Figure 2b.pdf", width=6, height=4)

rd.plot(rd.data$ever_winhouse, ever_winhouse)
ggsave("7tex/manuscript/tables/sourcefiles/Figure4a.pdf", width=6, height=4)

rd.plot(rd.data$ever_winsenate, ever_winsenate)
ggsave("7tex/manuscript/tables/sourcefiles/Appendix Figure 2d.pdf", width=6, height=4)

#win house
# rd.data %>% ungroup %>% 
#   filter(rd.data$victory_marg > -rd.data$bandwidth[1] & rd.data$victory_marg < rd.data$bandwidth[1]) %>% 
#   mutate(win = ifelse(victory_marg > 0, "win", "lose")) %>% 
#   ggplot(aes(x = victory_marg, y = ever_winhouse, group = win)) +
#   geom_smooth() +
#   geom_point(data = rd.data %>% 
#                filter(rd.data$victory_marg > -rd.data$bandwidth[1] & rd.data$victory_marg < rd.data$bandwidth[1]) %>% 
#                mutate(round.marg = round(victory_marg, digits = 2)) %>% 
#                group_by(round.marg) %>% 
#                summarize(hle.prob = mean(ever_winhouse)),
#              aes(x = round.marg, y = hle.prob), inherit.aes = F, alpha=.5) +
#   geom_vline(xintercept = 0, linetype = 2, color = "#666666", alpha=.8) + 
#   theme_bw() +
#   theme(panel.grid.minor.y=element_blank(),
#         panel.grid.minor.x = element_blank()) +
#   xlab("Victory Margin") + ylab("Prob. of Winning House Election")

  # geom_point(data = rd.data %>% 
  #              mutate(ntile = ntile(victory_marg, 30)) %>% 
  #              group_by(ntile) %>% 
  #              summarize(hle.prob = mean(ever_winhle)) %>% 
  #              mutate(victory_marg = seq(-.27, .27, length.out = 30)),
  #            aes(x = victory_marg, y = hle.prob), inherit.aes = F)



